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ABSTRACT 

We develop a model of the solar dynamo in which, on the one hand, we follow 
the Babcock-Leighton approach to include surface processes like the production 
of poloidal field from the decay of active regions, and, on the other hand, we 
attempt to develop a mean field theory that can be studied in quantitative detail. 
One of the main challenges in developing such models is to treat the buoyant rise 
of toroidal field and the production of poloidal field from it near the surface. The 
previous paper by Choudhuri, Schiissler, & Dikpati (1995) did not incorporate 
buoyancy. We extend this model by two contrasting methods. In one method, we 
incorporate the generation of the poloidal field near the solar surface by Durney's 
procedure of double ring eruption. In the second method, the poloidal field 
generation is treated by a positive a-effect concentrated near the solar surface, 
coupled with an algorithm for handling buoyancy. The two methods are found 
to give qualitatively similar results. 

Subject headings: MHD — Sun: interior — Sun: magnetic fields 

1. Introduction 

Historically there have been two theoretical approaches in understanding the origin 
of the solar magnetic cycle: the Parker-Steenbeck-Krause-Radler (PSKR) approach 
(Parker 1955; Steenbeck, Krause & Radler 1966) and the Babcock-Leighton (BL) approach 
(Babcock 1961; Leighton 1969). In both these approaches, the toroidal component of the 
magnetic field is supposed to be generated from the poloidal component by the stretching 
of field lines due to differential rotation. In order for a self-sustaining dynamo to exist, 
the poloidal field has to be generated back from the toroidal field. The two approaches 
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mentioned above attempt to solve this problem in two different ways. In the PSKR 
approach, the cyclonic turbulence in the interior of the Sun is supposed to twist the toroidal 
field lines to produce the poloidal field (the so-called a-effect). On the other hand, the BL 
approach puts more stress on what is happening at the solar surface and assumes that the 
poloidal field arises out of the decay of tilted bipolar active regions (which result from the 
magnetic buoyancy of the toroidal field). Various aspects of the generation of poloidal field 
at the surface have been discussed by Wang & Sheeley (1991) and Durney, De Young, & 
Roxburgh (1993). 

A formal mathematical formulation of the PSKR approach was developed on the basis 
of the mean field magnetohydrodynamics (Steenbeck, Krause, & Radler 1966; Moffatt 1978, 
Chap. 7; Parker 1979, §18.3; Choudhuri 1998, §16.5). In comparison, the BL approach 
was based on rather heuristic, and often qualitative, arguments. Until recently, most of 
the detailed mathematical models of the solar dynamo were worked out on the basis of the 
PSKR approach. Only in the last few years there have at last been attempts of putting the 
BL approach on a mathematical footing comparable in sophistication to the mathematical 
theory of the PSKR approach (Choudhuri et al. 1995; Durney 1995, 1996, 1997; Dikpati & 
Charbonneau 1999). It now appears that the most successful model of the solar cycle will 
be something which incorporates the best features of both these approaches (Choudhuri 
1999). 

Since magnetic buoyancy would be particularly destabilizing in the main body of the 
convection zone (Parker 1975; Moreno- Insert is 1983), several theorists (Spiegel & Weiss 
1980; van Ballegooijen 1982; DeLuca & Gilman 1986; Choudhuri 1990) argued that the 
solar dynamo may be operating in the overshoot layer at the bottom of the convection zone. 
With the helioseismic discovery of a shear layer at the bottom of the convection zone, it 
now appears fairly certain that the generation of the strong toroidal field by the stretching 
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of field lines must be taking place there. However, it seems unlikely that the whole dynamo 
process (as envisaged in the PSKR approach) occurs at the bottom of the convection zone. 
The studies of buoyant rise of the toroidal flux from there suggest that the toroidal field at 
the bottom of the convection zone should be of the order 10 5 G, substantially stronger than 
the equipartition value (Choudhuri & Gilman 1987; Choudhuri 1989; D'Silva & Choudhuri 
1993; Fan, Fisher, & DeLuca 1993; Caligari, Moreno-Insertis, & Schiissler 1995). Such a 
strong field would completely quench the a-effect of the PSKR approach. To explain the 
generation of the poloidal field, the most natural way is to invoke the BL idea of the decay 
of tilted active regions, though there are still some attempts to work within the PSKR 
approach by considering an interface dynamo (Parker 1993; Charbonneau & MacGregor 
1997; Markiel & Thomas 1999). In this paper, we assume that the poloidal field is produced 
by the decay of tilted active regions near the solar surface. 

Although the sunspots migrate equatorward with the solar cycle, the weak diffuse 
magnetic field on the solar surface migrates poleward (Bumba & Howard 1965; Howard & 
LaBonte 1981; Makarov, Fatianov, & Sivaraman 1983; Makarov & Sivaraman 1989). Most 
of the dynamo models based on the PSKR approach (starting from Steenbeck & Krause 
1969) mainly concentrated on the sunspots and ignored the poleward migration of the weak 
diffuse field. The poleward migration has been explained by assuming that the weak diffuse 
field (which is essentially the poloidal field) is carried by the meridional circulation (Wang, 
Nash, & Sheeley 1989a, 1989b; Dikpati & Choudhuri 1994, 1995; Choudhuri & Dikpati 
1999). If we now accept the BL idea that the poloidal field is produced by the decay of 
tilted bipolar active regions, then the meridional circulation should play an important role 
in the dynamo problem by bringing the poloidal field from the surface to the bottom of 
the convection zone, where the poloidal field is stretched out to produce the toroidal field. 
The challenge before us now is to develop a new type of dynamo model, in which the 
surface processes like the production of the poloidal field from the decay of active regions 
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are important as in the BL approach, but which has the same mathematical sophistication 
as the PSKR approach. Such a dynamo model should presumably account for both the 
equatorward migration of sunspots and the poleward migration of the weak diffuse field. An 
early step in this direction was taken by Wang, Sheeley, & Nash (1991), who averaged over 
the radial direction to obtain one-dimensional equations. More realistic two-dimensional 
models have been developed by Choudhuri, Schiissler, & Dikpati (1995), Durney (1995, 
1996, 1997) and Dikpati & Charbonneau (1999). 

The mathematical theory of the PSKR approach is based on mean field MHD, which 
leads to closed equations in the first order smoothing approximation. It is not clear if 
the implicit assumptions in this mathematical theory are fully satisfied in any realistic 
astrophysical situation. However, if the assumptions are satisfied, then the mathematical 
theory provides a completely rigorous description of the dynamo process in the PSKR 
approach. To make a similar rigorous formulation of the BL approach, we need to develop 
a consistent mean field description of (i) the buoyant rise of the toroidal flux to produce 
active regions and (ii) the decay of the tilted active regions to produce the poloidal field. In 
this paper, we focus our attention on comparing two possible formulations of the production 
of poloidal field from the decay of tilted active regions. Since it is necessary to include 
magnetic buoyancy to study this problem, we present some discussion of magnetic buoyancy 
as well. 

It was pointed out by Stix (1974) that the mathematical formulation of the BL 
approach is in some ways analogous to the a-effect of the PSKR approach. Choudhuri, 
Schiissler, & Dikpati (1995) modeled the decay of titled active regions to produce the 
poloidal field by invoking an a-coefficient which is concentrated near the solar surface. 
Durney (1995, 1996, 1997) followed Leighton (1969) more closely and treated the same by 
introducing a double ring of flux at the surface where the eruption takes place. Introducing 
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an a-coefficient concentrated near the surface is certainly a very approximate way of 
incorporating the main idea of the BL approach into the mathematical theory of the PSKR 
approach. Justifying this procedure rigorously is even more difficult than justifying the 
a-coefficient in the PSKR approach. However, this procedure produces the desired effect 
of generating the poloidal field where we want to generate it. If magnetic buoyancy is 
included in some way to bring the strong toroidal field from the bottom to the top and then 
the concentrated a-effect acts on it, the net result is similar to what happens in Durney's 
double ring method. Since this procedure is easier to implement than Durney's double ring 
method, one important question is whether this procedure is at least as good as Durney's 
double ring method. In this paper, we take a simple dynamo model and present calculations 
done with both the methods. We show that the results are qualitatively similar. It may be 
noted that it is not our aim to build realistic models of the solar cycle in this paper. For 
example, we have presented a contrasting study of these methods by assuming a differential 
rotation which does not vary with latitude as in Choudhuri, Schiissler, & Dikpati (1995). 
This simplification allows the specific features of the two methods to be seen clearly. A 
realistic differential rotation makes the results immensely more complicated, which we 
shall discuss in our next paper in which an attempt will be made to model the solar cycle 
properly. 

Durney (1995, 1995, 1997) allowed flux eruption to take place only at one latitude 
at a time. In Durney's model, it is difficult to allow simultaneous eruptions in a band of 
latitudes, which happens in the real Sun. The model of Choudhuri, Schiissler, & Dikpati 
(1995) did not incorporate magnetic buoyancy and allowed the toroidal field to be brought 
to the surface from the bottom by meridional circulation. To make comparisons with 
Durney's double ring method, we now include magnetic buoyancy in that model by allowing 
the magnetic field to erupt whenever it has a value larger than a critical value. It may be 
noted that incorporating magnetic buoyancy in the PSKR approach was relatively easier, 
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since magnetic buoyancy there merely removed the flux from the dynamo region and played 
the role of a dissipative process. Some authors treated magnetic buoyancy by putting a 
simple loss term in the dynamo equation (DeLuca & Gilman 1986; Schmitt & Schiissler 
1989), whereas others included a general upward flow due to magnetic buoyancy (Moss, 
Tuominen, & Brandenburg 1990a, 1990b). We have to go beyond such simple prescriptions 
in a BL approach, where magnetic buoyancy is a more integral part of the dynamo process 
and is not just a flux removal mechanism. In our BL model, magnetic buoyancy removes 
the flux from the bottom layer where the toroidal field is generated and then brings the 
flux to the top of the convection zone where the poloidal field is produced from it. Earlier, 
Choudhuri & Dikpati (1999) and Dikpati & Charbonneau (1999) incorporated the effect of 
magnetic buoyancy by including a dynamo source term near the surface which is a product 
of the a-coefficient and the toroidal magnetic field at the bottom of the convection zone. 

From the observation that the following spots in active regions appear at higher 
latitudes on the solar surface, it is easy to figure out that a has to be positive in the 
northern hemisphere. This is also clear from the expression of a-coefficient obtained by 
Stix [1974; Eq. (8)] by recasting the equations of Leighton (1969). The positive sign of a 
gives a new twist to the problem. It is well known that the product of a and the vertical 
gradient of differential rotation has to be negative in the northern hemisphere for the 
equatorward propagation of the dynamo wave (Parker 1955; see Choudhuri 1998, §16.6). 
Even if a and the velocity gradient are concentrated in two different layers, this condition 
still remains valid (Moffatt 1978, §9.7). Since the vertical gradient of differential rotation in 
the lower latitudes, as found by helioseismology, is positive, its product with a is positive 
and one would expect a poleward propagation of the dynamo wave. It was demonstrated by 
Choudhuri, Schiissler, & Dikpati (1995) that an equatorward propagation is still possible 
in this situation, if the time scale of meridional circulation is shorter than the time scale 
of diffusion between the layers of a and velocity shear. This opens up the possibility of 



building models of the solar dynamo in which we have a positive a near the surface and 
a positive gradient of differential rotation at the bottom of the convection zone. The 
meridional circulation has to play a very crucial role in such models in ensuring the desired 
behaviour. While using Durney's double ring method, the signs of the magnetic field in 
the two rings have to be chosen such that there is a correspondence with the positive a 
situation. With the double ring method also, we found that the dynamo wave at the bottom 
of the convection zone propagates equatorward only when there is a strong meridional flow 
and propagates poleward when this flow is switched off. 

In §2 we discuss the details of our model. Then we go on to present our main results 
in §3. Our conclusions are summarised in §4. 

2. The model 

We assume axisymmetry in all our calculations. The magnetic and velocity fields can 
be written as 

B = Be 4> + Vx{Ae 4> ), (1) 

v = v p + r sin OVte^. (2) 

where B and A respectively represent the toroidal and poloidal components of the magnetic 
field; Q is the angular velocity, and v p = v r e r + v geo is the meridional circulation. We 
substitute equations (1) and (2) in the induction equation 

— = Vx(vxB)+7?V 2 B, (3) 

where 77 is the coefficient of turbulent diffusion. This gives 

|V(v,V)M)=r?(V-l).4 + Q, (4) 



dB 1 
dt r 



|(nvfl) + |(«^)' 



>l (V 2 -^)B + s(B,,.V)i» iV) 
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where s = rsin# and B p = V x (Ae^). We have added one extra term Q on the right-hand 
side of equation (4), which does not follow from the induction equation (3). It is a term 
which describes the generation of the poloidal field. The usual afi dynamo is given by the 
equations (4) and (5), where Q is simply 

Q = aB. (6) 

To incorporate the effect of magnetic buoyancy and the decay of tilted active regions, 
we have to allow for changes in B due to the rise of magnetic flux from the bottom of 
the convection zone to the top and specify Q appropriately. Before describing how we 
incorporate Durney's double ring method as well as our method of concentrated a-effect 
near the surface, let us discuss a few general points which hold for both cases. 

The equations (4) and (5) have to be solved in the northern quadrant of the convection 
zone as usual (i.e. within Rf, = 0.7R Q < r < R Q , < < 7r/2). The boundary conditions 
are discussed in previous papers (Dikpati & Choudhuri 1994; Choudhuri, Schiissler, & 
Dikpati 1995). They are 

At 6 = : A = 0, B = 0, 

d 
Atr = R b : A = 0, -z-irB) = 0, 

or 

Atr = R Q : B = 0, 

the boundary condition for A at the top r = R & being that it has to match a smooth 
potential field outside. See §3 of Dikpati & Choudhuri (1994) for a detailed discussion of 
how this is implemented. 

To solve equations (4) and (5) with these boundary conditions, we need to specify r), fl, 
v p and Q. As in Choudhuri, Schiissler, & Dikpati (1995), we assume the turbulent diffusion 
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to have the constant value 77 = 1.1 x 10 7 m 2 s _1 . For the angular velocity Q also, we use 

the same expression as used in that paper: 

- r — r3 



1 + erf 



rfa 



} 



(7) 



fi = fi {0.9294 + 0.0353 

with r 3 = 0.7-Rq, d 3 = 0.1R Q , Q = 2.7 x 10~ 6 s _1 . This latitude-independent angular 
velocity roughly corresponds to the helioseismologically determined rotation profile near 
the solar equator, with dVt/dr positive. For the meridional circulation v p , we again use the 
expression used previously (Dikpati & Choudhuri 1995; Choudhuri, Schiissler, & Dikpati 
1995). In other words, we take 

pv p = V x (>e ) (8) 



n{r - R b 



{1 _ e _/3ir9E } x {1 - e /32r(e-7r/2) j e -((r-n,)/r) 2 



with ip given by 

ipr sin 9 = ipQ sin 
and p given by 

The values of the parameters used are f3\ = 1.4 x 10~ 8 m _1 , fa = 2.7 x 10" 



p(r) = C[^-~ 



(9) 



(10) 



m 



-1 



e = 2.0000001, r = (R Q - R b )/5, T = 3.3 x 10 8 m, 7 = 0.9, m = 3/2. The pattern 
of meridional circulation for these values of parameters is shown in Figure 3a of Dikpati 
& Choudhuri (1995). The amplitude of the meridional circulation is fixed by taking 
ipo/C = —7.9 x 10 8 m 2 s _1 , which corresponds to a maximum surface velocity (i>o) of about 
7.0 m s _1 in the mid-latitudes. 

In the paper of Choudhuri, Schiissler, & Dikpati (1995), the source term Q was given 
by equation (6) with a taken in the form 



a 



«o 



5 2COS 4 



1 + erf 



r — r\ 
d\ 



x 



erf 



r — T2 
do 



'IF 



The parameters are r\ = O.95i? , r2 = Rq, d\ = di = 0.025R & , making sure that the a 
effect is concentrated in the top layer O.95i? < r < Rq. The a-quenching factor 1 + B 2 
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included in the denominator helps the system to relax to periodic solutions with amplitude 
B ~ 1. This essentially means that we are choosing the unit of B in such a way that the 
nonlinear feedback becomes important when B is of order unity or larger. Choudhuri, 
Schiissler, & Dikpati (1995) took Qo = 3m s _1 and found that it gave rise to marginally 
critical oscillations. When magnetic buoyancy is included, we find that this value of «o 
often gives decaying solutions. To ensure that the solutions do not decay, we take a = 10 
m s _1 in most of the calculations in the present paper. 

2.1. Incorporating the double ring 

After time intervals r, we find the co-latitude 8 er where the toroidal field is maximum 
and allow the flux to erupt above in the form of the double ring, if this maximum value 
exceeds a specified critical field B c . Following Figure 1 of Durney (1997), we show the 
two emergent flux rings in Figure 1. One ring of positive magnetic field K/ sm.8 is put 
between the co-latitudes 0±, 8 2 , whereas the other ring of negative magnetic field — K/ sind 
is between #3, 84. The factor sin 8 ensures that the flux through one ring balances the flux 
through the other ring. In Durney's notation, 8±, 82, #3 and #4 will be 

8\ = 8 CI — 

82 = 8 cr — 

83 = 8 CT + 

84 = 8 cr + ; 

As in Durney (1997), we make the somewhat unphysical assumption that these rings extend 
only from R Q to R & — Ar, where the field lines end abruptly. At the time of eruption, then, 
in the region R Q — Ar < r < R & , we put the magnetic field 
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AB r = if 9i < 9 < 9 2 

sin 9 

~ h if 9, <9<9 d 



sin 9 
= elsewhere (12) 

Putting this magnetic field is equivalent to adding the vector potential AA given by 

AB r = —^- — (sin 9 A A), 
r sin 9 39 

from which 

r r e 

AA = -/ sin9'AB r d9' (13) 

sin 9 J 

if we do not consider the variation of AA in r. Substituting for AB r from (12), we conclude 
that AA can be non-zero only in the range R Q — Ar < r < R Q , where we have 

AA sin 9 =0 for < 9 < 9 X 

= R Q K(9 - 0x) for 9 ± <9<9 2 

= R Q K(9 2 - 0±) for 9 2 <9<9 3 

= R Q K[(9 2 - Ox) -{9- 9 3 )] for 9 3 < 9 < 9 A 

= for 9 4 < 9 < - (14) 

Adding this AA to A leads to a discontinuity in A at R Q — Ar. Durney (1997) writes, "Such 
an expression for the vector potential generates latitudinal magnetic fields (associated with 
the closure of magnetic lines of force)", but also claims that these discontinuities "are 
numerically inconsequential" . 

Durney (1995, 1996, 1997) took the separation between the rings, \i to be proportional 
to cosher- To keep the numerical computations simpler, we instead take K appearing in 
(14) to be proportional to cos# er , while keeping the separation between the rings fixed. 
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This has the same physical effect, except at the low latitudes where the two rings may 
overlap. We, however, find that flux eruption remains restricted to higher latitudes where 
this overlap is unimportant. If we take K to be proportional to the toroidal magnetic field 
B at the bottom (from which the flux rings originate), then our problem becomes linear in 
magnetic field and one has to make many runs to find the marginally growing solution. We 
circumvent this problem by including something like a-quenching in the following fashion: 

K = K i + \ B r (15) 

■J- T^ | -'-''max | 

where -B max is the toroidal magnetic field at the bottom of the convection zone at the 
latitude where it is maximum. The justification behind this is the fact that a stronger 
toroidal field is less affected by the Coriolis force (D'Silva and Choudhuri 1993; Howard 
1993) and hence is less efficient in generating the poloidal field. It can easily be seen from 
(14) and (15) that K' is a dimensionless quantity. 

It is seen on the solar surface that the magnetic field of the higher-latitude sunspot 
is positive when the toroidal field underneath the surface is positive. It should be clear 
from (12) and (15) that this is achieved by taking K' positive, which is the case in all our 
calculations. It now follows from (14) and (15) that a positive B at the bottom of the 
convection zone would imply a positive increment in A at the surface where the magnetic 
flux emerges. This is something like a positive a-effect, which obviously corresponds to a 
positive value of K' . 

We solve (4) and (5) with differential rotation and meridional circulation as given by 
(7)-(10). The source term Q in (4) is given by (6) and (11). However, in addition to this 
usual source term, we allow for possible changes in the value of A abruptly, in the double 
ring regions of the surface, at intervals of r, to take account of magnetic buoyancy. We run 
our code to find the maximum value of B after intervals r. If this exceeds the critical field 
B c and occurs at the co-latitude 8 er , then we consider two rings situated on two sides of this 
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co-latitude and add AA as given by (3) to A. The control parameter in our problem is K' 
appearing in (4). When K' is zero, there is no double ring formation and we get the model 
of Choudhuri, Schiissler, & Dikpati (1995). On the other hand, when K' is sufficiently large, 
the net effect of double ring formation at intervals of r becomes much more important 
than the source term Q in (4) and we have the model of Durney (1997). Thus, in the two 
opposite limits of the control parameter K', our model is respectively reduced to the models 
of Choudhuri, Schiissler, & Dikpati (1995) or Durney (1997). 

2.2. Concentrated a-effect with magnetic buoyancy 

We wish to argue that the double ring method is similar to allowing magnetic flux to 
rise due to magnetic buoyancy and then letting the a-effect concentrated near the surface to 
act on it. In this method also, we solve (4) and (5) in conjunction with (6)— (11). However, 
instead of having double ring formations at intervals of r (leading to abrupt changes of A 
as seen from [14]), we now allow B to change abruptly at intervals of r to take account of 
flux rise due to magnetic buoyancy. This is done in the following way. 

We assume that the toroidal field B becomes buoyant when its value crosses a critical 
value B c . After intervals of time r, we check if B has become larger than B c at certain 
points. Then, at those points, B is reduced by a factor 1 — /, i.e. 

fl-£(l-/). 

The flux removed from these points is taken vertically above and deposited near the 
surface by increasing B there in such a fashion that the total flux remains conserved in the 
transfer process. Since the equations are numerically solved on a JV x M grid, the simplest 
procedure is to deposit all the flux at the grid point just below the surface. For example, 
if B crosses B c only at one grid point on the radial line at a fixed latitude, then we have 
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to decrease B by fB there and the toroidal field at the grid point just below the surface 
has to be increased by an amount f'B. Since the grid size at the surface corresponds to 
a greater distance in the latitudinal direction than that at the bottom, we need to take 
/' = f(Ri/Rf) to ensure the conservation of magnetic flux (here, Ri is the radius near the 
bottom where the flux is depleted and Rf is the radius near the surface where the flux is 
deposited). We have also made some runs in which the flux taken up from one grid point 
is distributed within a few grid points near the top instead of all the flux being deposited 
in one grid point, and the results turn out to be qualitatively similar. The strength of 
magnetic buoyancy is increased by increasing the control parameter /. In the limit / = 0, 
we get back the model of Choudhuri, Schiissler, & Dikpati (1995), in which there was no 
magnetic buoyancy and the toroidal field was brought to the surface by the meridional 
circulation. When / is made sufficiently large (even though it has to remain less than 1), 
magnetic buoyancy is found to dominate and the system has a limiting behavior. 

Compared to the double ring method, this method has some attractive features. 
Firstly, here the eruption at any instant takes place over a range of latitude rather than 
at one point as in the double ring method. This corresponds to the real Sun more closely. 
It is not easy to extend the double ring method to handle simultaneous flux eruptions at 
more than one point. If we simultaneously put several double rings in a range of latitudes, 
then the positive ring of an intermediate double ring will cancel with the negative ring 
of the next double ring and we shall be left with a positive ring and a negative ring at a 
wide separation. It follows from (14) that this will mean adding to A over a wide range of 
latitude. This would make the model more similar to the mean field model and the special 
character of the original double ring model would be completely lost. Also, we now allow 
for the toroidal flux to be depleted at the bottom of the convection zone due to magnetic 
buoyancy. As we shall argue later, we believe this to be quite important. In fact, we shall 
present some results with the double ring method with the toroidal flux at the bottom 
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depleted parametrically. 

3. Results 

We now present and compare results obtained by the two methods described above. As 
we saw, K' and / happen to be the respective control parameters in these two methods. On 
setting these control parameters equal to 0, both these methods are reduced to the model 
of Choudhuri, Schussler, & Dikpati (1995, hereafter CSD model). All our calculations are 
done on a 64 x 64 grid. We allow the eruptions to take place after times r = 8.8 x 10 5 
s and use a value B c = 1 for the critical field in all our calculations. When we start our 
calculations with any arbitrary magnetic field configuration, the code relaxes to a periodic 
solution for a proper set of parameters. What we discuss below are properties of such 
relaxed periodic solutions. 

3.1. Results with the double ring method 

Durney (1997) did not allow the toroidal flux to be depleted at the bottom of the 
convection zone due to magnetic buoyancy. To study the effect of flux depletion, we present 
some calculations in which we allow flux depletion in the following simple manner. At the 
times of eruption after interval r, we find out at which point the toroidal field has the 
maximum value B max (> B c ). While putting the two flux rings at the top, we also decrease 
-B max by an amount fd,B max at the maximum point. Then fd becomes a second parameter 
in the problem in addition to K 1 in our problem. After finding the co-latitude 8 CT where the 
toroidal field is maximum, the next two poleward grid points are taken as $i, $ 2 , and the 
next two equatorward grid points are taken as #3, #4. The flux rings are assumed to go 3 
grid points deep (i.e. Ar is taken 3 grid points below the surface). 
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Figure 2 shows how the dynamo period Td changes with the parameter K' when fd is 
held constant. The different curves correspond to different values of fd- When we go to 
the limit of CSD model by putting K 1 — 0, we find the period to be 66 yrs. When fd — 
(i.e. there is no flux depletion at the bottom), we find that the change in the period with 
K' does not follow any particular trend. Td at first increases slightly with increasing K' 
and then comes down to a value close to that of the CSD model. This behaviour for fd — 
may result from the fact that in this case we are actually creating flux (in the form of 
erupted double rings) without any depletion. More meaningful behaviour follows for the 
other values of fd (such as 0.25, 0.5, 0.75). The period decreases with increasing K' and 
tends to saturate at some asymptotic value for large K' . To understand what is happening, 
let us look at Figure 3 which shows the evolution of magnetic field during a half-period 
for the case K' = 1000, fd = 0.5. In the plots of poloidal field, we have indicated the 
latitudes of last flux eruption with small arrows. However, the individual double rings are 
not usually discernable. That is not surprising. Flux eruption in the form of double rings 
keeps occuring at intervals of r. Hence the latest double ring is merely superposed on the 
field created by the previous double rings and does not stand out against the background 
of previously created field. On looking at the plots of the toroidal field, it is clear that 
the toroidal field keeps weakening as we go to lower latitudes. This weakening of toroidal 
field at lower latitudes becomes more prominent as we make fd larger. This implies that 
flux eruption never takes place at very low latitudes and the dynamo process is basically 
confined to higher latitudes. Since it takes less time to transport magnetic flux through a 
limited range of latitudes, the dynamo period is shorter for non-zero fd- In combination 
with this effect, an increasing K' will make the erupted double rings stronger, thus recycling 
toroidal flux to poloidal flux more efficiently. This reduces the time period of the dynamo 
as compared to the period in the limit of the CSD model, in which the toroidal field is 
brought to the surface by the meridional flow only near the equator and the whole range 
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of latitudes is involved. It may be noted that Durney (1997) did not present any plots of 
magnetic field configurations in his paper. However, we do get a deeper insight into the 
problem by looking at such field configuration plots. For example, note that the direction 
of the poloidal field (clockwise or anti-clockwise) starts reversing at the time when we have 
an extended belt of strong toroidal field. 

Durney (1997) has presented several plots showing how the eruption latitude changes 
with time (Figures 7-10 in his paper). We present a similar plot in Figure 4 for the case 
K' = 1000, fd = 0, corresponding to no flux depletion at the bottom as in the calculations 
of Durney (1997). Here, we see that eruptions continue near the pole for some time at 
the beginning of a cycle and then progressively move to lower latitudes. This plot looks 
very much like the plots presented by Durney (1997) -- especially his Figure 7. This is 
certainly very reassuring, since the numerical techniques employed by us and by Durney 
(1997) are completely different. Apart from the production of the double rings, our code 
allows for the toroidal flux to be brought to the surface by meridional circulation and then 
to be acted upon by a-coefficient (an effect not present in Durney's calculations). However, 
when K' is made as large as 1000, this effect is insignificant. In fact, we made some runs 
with a = and found that the results for zero or non-zero a are virtually indistinguishable 
when K' = 1000. For example, the plots of eruption latitude against time and the butterfly 
diagrams look identical in both the cases. 

We have already mentioned that a positive K' is like a positive a-effect concentrated 
near the surface. Choudhuri, Schussler, & Dikpati (1995) showed that a positive a 
concentrated near the surface leads to a poleward propagation of the dynamo wave when the 
meridional flow is switched off. We find exactly the same result in the double ring approach 
with positive K' if we switch off the meridional flow. Figure 5 shows a time-latitude plot 
of the toroidal field at the bottom of the convection zone with meridional flow for the case 
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K' = 1000, fd = 0.5, whereas Figure 6 is a similar plot without meridional flow keeping all 
the other parameters the same. We see clear indication of poleward migration in Figure 6. 

3.2. Results for concentrated a with buoyancy 

For contrast, we now present results obtained by the method described in § 2.2. As we 
have seen, the control parameter in this problem is /(< 1), which measures the strength of 
magnetic buoyancy. Figure 7 shows how the dynamo period changes on increasing /. As 
in Figure 2, we begin with a period of 66 yrs in the limit / = corresponding to the CSD 
model. On making the effect of buoyancy stronger (by increasing /), the flux transport 
(from the bottom of the convection zone to the top) takes place more efficiently and also 
the toroidal flux gets depleted quickly. This results in the dynamo period reducing with 
increasing /, until it reaches an asymptotic value of about 25 yrs. We may point out here 
that we did some runs for this method without depleting the field at the bottom, which 
would correspond to the case fd — for the double ring method. We found that even in 
this case, there is no decrease in the time period with increasing / (which in this case 
corresponds to only field addition at the top) and Td more or less hovers around the CSD 
limit of 66 years. 

Since the two methods discussed by us are sufficiently different, it is not obvious 
which value of K' in the first method would correspond to a certain value of / in the 
second method. In both methods, however, the dynamo periods saturate to asymptotic 
values when these control parameters are sufficiently large. So the most sensible thing is 
to compare results of the two methods when the control parameters are large enough to 
ensure that the dynamo period has the asymptotic value. Figure 8 shows the time evolution 
of the magnetic field during a half period for the parameters / = 0.05 (i.e. the magnetic 
buoyancy is strong enough to saturate the period to its asymptotic value). On comparing 



-20- 

with Figure 3, we find that the broad features of the magnetic field distribution are very 
similar. The main difference is that one sees some toroidal field distributed near the top of 
the convection zone in Figure 8, whereas such fields are not present in Figure 3. The reason 
behind this is obvious. In the double ring method, we directly put double rings above 
regions of strong toroidal field and this contributes directly to the poloidal field. When 
we introduce the intermediate step of the toroidal field first rising due to buoyancy and 
then being acted upon by a-effect, then we get toroidal field at the top of the convection 
zone also, as in Figure 8. The other difference between Figures 3 and 8 is that often the 
field lines in Figure 3 in some places (especially near the surface) are not as smooth as 
the field lines are everywhere in Figure 8. This is certainly due to double ring formations 
in Figure 3, which are concentrated local effects. As in Figure 3, here also we find that 
the direction of poloidal field reverses at around the time the strong toroidal field belt is 
maximally extended. 

Finally Figure 9 presents a time-latitude plot of the toroidal field at the bottom for the 
same case which is presented in Figure 8. Again, this figure looks qualitatively similar to 
Figure 5, the main difference being the fact that the toroidal field has become much weaker 
near the equator in Figure 8 due to more efficient flux depletion at the bottom, which takes 
place naturally in this method. 

4. Conclusion 

Following Choudhuri, Schiissler, and Dikpati (1995) and Durney (1997), we build a 
hybrid model of the solar dynamo, in which the best features of both the PSKR and the BL 
approaches are combined. The aim is to include the surface processes emphasized in BL 
models into a model as suitable for detailed quantitative study as the PSKR models. We 
study two possible methods of achieving this. One is to introduce double rings above the 
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region where the toroidal field is maximum, as done by Durney (1995, 1996, 1997). The 
second method is to make the toroidal field rise when it is above a critical value and then 
allow it to be acted upon by an a-coefficient concentrated near the surface. It is reassuring 
that the results obtained by the two methods are qualitatively similar. 

We believe that the depletion of toroidal flux by magnetic buoyancy is an important 
process. Flux tube calculations (Choudhuri & Gilman 1987; Choudhuri 1989; D'Silva 
& Choudhuri 1993; Fan, Fisher, & DeLuca 1993; Caligari, Moreno-Insertis, & Schiissler 
1995) suggest that the toroidal field at the bottom of the convection zone has a value of 
10 5 G — much stronger than the equipartition value. After the belt of strong toroidal field 
reaches the equator, it disappears and the next half-cycle of the dynamo begins. If the field 
is so strong, then turbulent diffusion will be completely suppressed and will not be effective 
in destroying the strong toroidal field. The only way to annihilate this belt of strong 
toroidal field is to expect magnetic buoyancy to deplete its strength sufficiently by the time 
this belt propagates to the equator. In our second method, this flux depletion automatically 
takes place. In the double ring method, we have included the possibility of toroidal flux 
depletion as an extra effect, which was not taken into account by Durney (1997). When the 
toroidal flux is depleted appropriately, both the methods make the period of the dynamo 
decrease on increasing the control parameters (K f or /) and saturate at some asymptotic 
value. This decrease of period is due to the efficient and rapid transport of toroidal flux 
by magnetic buoyancy. It is true that the decrease of period is more pronounced in the 
second method, where the flux depletion is more prominent. However, a difference by 
a factor 2 or 3 in the asymptotic period is probably not such a significant uncertainty 
compared to many other factors. The magnetic field configurations obtained in the two 
methods, as seen in Figures 3 and 8, are also quite similar, with the poloidal field reversing 
its direction for the same configuration of the toroidal field. Then, very importantly, the 
dynamo wave is found to propagate poleward when the meridional circulation is switched 
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off in the double ring method. In other words, the double ring method with positive K' has 
characteristics quite similar to a model with positive a-effect concentrated near the surface. 
The results obtained by the two methods are not exactly identical. However, given the 
many uncertainties plaguing the solar dynamo theory at the present time, representing the 
generation of poloidal field near the surface by a concentrated a-effect acting on erupted 
toroidal field seems like a good enough approximation. 

We should point out that there are several logistic problems in numerically handling 
the double ring, which are not there if we use an a-effect instead. Firstly, to properly 
create rings of latitudinal size similar to sunspot size with appropriate separation, one has 
to either use at least 500 grid points in the direction or use a special code which employs 
a finer mesh in the region where eruption takes place. Durney (1997) used 101 x 101 grid, 
which corresponds to a grid size of about 11,000 km in the latitudinal direction at the Sun's 
surface. The width of the double ring has to be at least 4 times this, i.e. about 44,000 
km — definitely inadequate to resolve the north-south polarity separation of a typical active 
region. To ensure whether our results have converged with respect to grid size, we repeated 
some calculations on 32 x 32 grid and compared the results with those obtained on our usual 
64 x 64 grid. We found that results obtained by our second method of concentrated a-effect 
were so close in the two cases that various plots looked indistinguishable. However, results 
obtained by the double ring method, in which important source terms are taken at the limit 
of grid resolution, while remaining qualitatively similar on halving the grid size, showed 
some changes. Our grid size is comparable to what other researchers (Durney 1997, Dikpati 
& Charbonneau 1999) have used on similar problems. We believe that the grid size has 
to be reduced considerably to properly resolve double rings and to give results completely 
invariant with grid size. Since the ring separation was at the limit of grid resolution, we 
kept the ring separation fixed and made our constant K proportional to cos# er (see [15] 
and the discussion preceding that). Durney (1997) claims to have made the ring separation 
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proportional to cos# cr , but never explains in his paper how this could be done with only 
101 grid points in the 9 direction. Another important consideration is that the double ring 
method is easy to implement when we allow flux eruption only at one point at one time, but 
it is not easy to generalize if multiple flux eruptions are allowed. In reality, we find that, 
at a certain time, several active regions emerge in a belt of latitudes — with the different 
active regions usually separated in longitude. If one could use an appropriately resolved 
3D code in which active regions of realistic size were made to emerge in different latitudes 
and longitudes, then certainly that would have been a much more satisfactory calculation 
than what we are doing now. We hope that future computers will be used by researchers 
more numerically capable than us to tackle this problem. When one uses rings to replace 
active regions through an averaging over longitude and uses a grid not fine enough to 
resolve individual sunspots, one already introduces some drastic averaging. Introducing an 
a-coefficient concentrated near the surface instead of using double rings may not be such a 
big step after that. 

Let us end with a comment on what we mean by a Babcock-Leighton model, since this 
term often creates some confusion. Babcock (1961) and Leighton (1969) emphasized the 
surface process of poloidal field generation from tilted active regions — in contrast to the 
usual mean field MHD where the poloidal field is supposed to be produced in the interior 
region of turbulence (see, for example, Choudhuri 1998). Hence any dynamo model in 
which the poloidal field is generated in a thin layer near the solar surface should be called a 
Babcock-Leighton model. Dikpati & Charbonneau (1999) also use the term in this sense. 
Durney (1995, 1996, 1997) followed Leighton (1969) more closely in incorporating the 
Babcock-Leighton idea through the double ring method. Introducing a phenomenological 
a-coefficient concentrated near the surface is another way of representing the Babcock- 
Leighton idea. One should, however, be careful not to interpret this a-coefficient in the way 
it is interpreted in the mean field MHD. For example, the a-coefficient here is not obviously 
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related to the average helicity of turbulence as in mean field MHD (see, for example, 
Choudhuri 1998, §16.5). This coefficient merely provides a phenomenological description of 
the production of poloidal field from the decay of tilted active regions, which is obvious in 
the formulation of the BL model by Stix (1974). Wang & Sheeley (1991) also referred to 
this process as an "a-effect" in exactly the same sense as us, even though they never used 
the symbol a in their actual equations! 

We would like to thank Paul Charbonneau, Bernard Durney, Gene Parker and an 
anonymous referee for valuable suggestions. 
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Fig. 1. — A sketch of the erupted double ring at the solar surface. 
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Fig. 2. — The variation in the dynamo period (in units of years) with the control parameter 
K f for four different fd values. The dash-dotted line corresponds to fd — 0, the solid line to 
fd = 0.25, the dotted line to fd = 0.5 and the dashed line to fd = 0.75. 
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Fig. 3. — Time evolution of the toroidal field (left hand column) and poloidal field (right hand 
column) configuration in a meridional cut of the northern quadrant of the solar convection 
zone (O.7i? , < r < R e , < 9 < vr/2) for the case with K' = 1000, f d = 0.5. The whole set 
covers a dvnamo half period. That is from toD to bottom t = 0. Tw/8. Tw/4. 3T^/8. TaII. 
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Fig. 4. — The variation in the eruption latitude with time (in years) for K' = 1000, fd = 0. 




Fig. 5. — Time-latitude plot of the contours of constant toroidal field B at the bottom of 
the convection zone, for the case with K' = 1000, fd = 0.5 and vq = 7.Qms~ l . The solid 
lines denote positive B and the dashed lines negative B. Time is in years. 
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Fig. 6. — Time-latitude plot of the contours of constant toroidal field B at the bottom of 
the convection zone, for the case with vq = 0.0. The other parameter are the same as in 
Figure 5. 
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Fig. 7. — The variation in the dynamo period (in units of years) with the control parameter 
/ for our second method - concentrated a effect with buoyancy. 
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Fig. 8. — Time evolution of the toroidal field (left hand column) and poloidal field (right hand 
column) configuration for the concentrated a effect with buoyancy method with / = 0.05. 
The convention followed is the same as in Figure 3. 
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Fig. 9. — Time-latitude plot of the contours of constant toroidal field B at the bottom of the 
convection zone for the concentrated a effect with buoyancy method, with / = 0.05. Time 
is in years. 



